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A 'forward walking' Green's Function Monte Carlo algorithm is used to obtain 
expectation values for SU(3) lattice Yang-Mills theory in (3+1) dimensions. The 
ground state energy and Wilson loops are calculated, and the finite-size scaling 
behaviour is explored. Crude estimates of the string tension are derived, which 
agree with previous results at intermediate couplings; but more accurate results 
for larger loops will be required to establish scaling behaviour at weak coupling. 



1 Introduction 

Classical Monte Carlo simulations provide a very powerful and accurate method 
for the study of Euclidean lattice gauge theories. In the Hamiltonian formula- 
tion, on the other hand, the corresponding quantum Monte Carlo methods have 
been somewhat neglected. Here we present a study of SU(3) Yang-Mills theory 
in (3-1-1) dimensions, using the Green's Function Monte Carlo approachEl 

Heys and StumpB and Chin et al. □ pioneered the use of "Green's Function 
Monte Carlo" (GFMC) or "Diffusion Monte Carlo" techniques in Hamiltonian 
LGT, in conjunction with a weak-coupling representation involving continuous 
gauge field lirLkpyariables. This was successfully adapted to non-Abelian Yang- 
Mills theoriesp'u with no minus sign problem arising. In this representation, 
however, one is simulating the wave function in gauge field configuration space 
by a discrete ensemble or density of random walkers: it is not possible to 
determine the derivatives of the gauge fields for each configuration, or to enforce 
Gauss's law explicitly, the ensemble always relaxes back to the ground state 
sector. In order to compute the string tension or mass gap, one must measure 
an appropriate correlation function, and estimate the mass gap as the inverse 
of the correlation length. We haveintroduced the 'for ward- walking' technique, 
well-known in many-body theoryp'tl to measure the expectation values and 
correlation functions. The technique Jias been demonstrated for thejj^ases of 
the transverse Ising model in (1+1)00 and the U(l) LGT in (2+1)dI 

Here we apply the technique for the first time to a non-Abelian model, 
namely SU(3) Yang-Mills theory in (3+l)D. The grotmd state energy and 
Wilson loop values are calculated, and approximate values are extracted for 
the string tension in the weak-coupling regime. Comparisons are made with 
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earlier calculations, where they are availablel 



2 Method 

2.1 Lattice Hamiltonian 

The Green's Function Monte Carlo formalism has been adapted to SU(2) Yang- 
Mills theory by Chin, van Roosmalen, Umland and Kooninp and sketched for 
the SU(3) case by Chin, Long and RobsonB ^ 
The SU(3) lattice Hamiltonian is given byQ 

I V 

where Ef is a component of the electric field at link 1, A = 6/5^, the index a 
runs over the 8 generators of SU(3), and Up denotes the product of four link op- 
erators around an elementary plaquette. We will work with the dimensionless 
operator 

H^\Y. EtEf -^Y. + (2) 

I p 

The link variables are elements of the group SU(3) in the fundamental 
representation 

U = exp(-4A"A'') (3) 



2.2 Green's Function Monte Carlo method 

The Green's Function Monte Carlo method employs the operator exp(— T(iJ — 
E)), i.e. the time evolution operator in imaginary time, as a projector onto the 
ground state IV'o): 

IV'o) oc lim e-^^"(^-^)|$) (4) 

Ar— >0,Af At-^oo 

where |<i>) is any suitable trial state. To procure some variational guidance, 
one performs a "similarity transformation" with the trial wave function $, and 
evolves the product $|'(/'o) i^i imaginary time. The heart of the procedure is the 
calculation of the matrix element corrersponding to a single small time step 
Ar. Chin et ala show that 

(x'|$e-^-(«-^)a>-i|x)=n(^/l^{^^P(-^^^^r£^r)exp[Ar£;r(£;rin<i>)]}|C/z) 

exp{AT[£; - $"ii/$(x)]} + 0(At2) 
= p(x',x)w(x) -|-0(Ar2) (5) 
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where x — {Ui} denotes an entire lattice configuration of link fields. 

The product ^\ip) is simulated by the density of an ensemble of random 
walkers. At the kth. step, the 'weight' of each walker at is multiplied by 
w{'}Ck). The effect of p(xk+i,Xk) is to alter each link variable Ui in {xk} to 
U[ by a Gaussian random walk plus a "drift step" guided by the trial wave 
function: 

U' = AUUdU (6) 

where Ud = exp[iiA"(iAT_E'' ln<i>)] is the drift step, and AU is an SU(3) group 
element randomly chosen from a Gaussian distribution around the identity, 
with variance (As^) = 8At (i.e. At for each index a), where 

(As2) « = 8Ar, (7) 

a 

for small A". 

The simulation is carried out for a large number of iterations At, until 
an equilibrium distribution <I>|V'o) is reached. The energy E in (^) is adjusted 
after each iteration so as to maintain the total ensemble weight constant. The 
average value of E can then be taken as an estimate of Eq, the ground-state 
energy. 

As time evolves, the weights of some walkers grow larger, while others 
grow smaller, which would produce an increased statistical error. To avoid 
this, a "branching" process is employed, whereby a walker with weight larger 
than some threshold is split into two independent walkers, while others with 
weights lower than another threshold are amalgamated. 

2.S Trial Wave Function 

The trial wave function is chosen to be the one-parameter formi 

$ = exp[a^rr(C/p + C/;)] (8) 
p 

Then the drift step for each linkisB 

C/rf = exp[-z^^f], (9) 

Af = -iAr^ J2 TriX^'Ui-Ul - h.c] (10) 
pel 
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Finally, the trial energy factor is 

$-ii/a> = ^{^(^Tr[A'^[/,..[/l - h.c.]f (11) 

I pel 

+(^-^-^)Y.T^^Up + ul)"^- (12) 

pel 

2.4 Forward Walking estimates 

The "forward walking" technique is used to estimate expectation valuesB Its 
application to the U(l) lattice gauge theory in (2+l)D was discussed by Hamer 
et alQ It is implemented for an operator Q (assumed diagonal, for simplicity) 
bycl recording the value (5(xi) for each "ancestor" walker at the beginning of a 
measurement; propagating the ensemble as normal for J iterations, keeping a 
record of the "ancestor" of each walker in the current population; and taking 
the weighted average of the Q(xi) with respect to the weights of the descen- 
dants of Xi after the J iterations, using sufficient iterations J that the estimate 
reaches a 'plateau'. 

3 Results 

Simulations were carried out for LxLxL lattices up to L=8 sites, using runs 
of typically 4000 iterations and an ensemble size of 250 to 1000 depending 
(inversely) on lattice size. Time steps At of 0.01 and 0.05 "seconds" were used, 
with each iteration consisting of 5 sweeps and 1 sweep through the lattice, 
respectively, followed by a branching process. The first 400 iterations were 
discarded to allow for equilibration. The data were block averaged over blocks 
of up to 256 iterations, to minimize the effect of correlations on the error 
estimates. 

The results taken at Ar = 0.01 and At — 0.05 were extrapolated linearly 
to A|i = 0. The variational parameter c was given values as used by Chin 
et aip obtained from a variational Monte Carlo calculation. We checked that 
these were approximately the optimum values for small lattices. 

Forward-walking measurements were taken over J iterations, where J 
ranged from 20 to 100, depending on the coupling A. Ten separate measure- 
ments were taken over this time interval, in order to check whether the value 
measured by forward- walking had reached equilibrium. A new measurement 
was started soon after the previous one had finished. 
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3.1 Ground-state Energy 

The dependence of the ground-state energy per site on lattice size is ilhistrated 
in Figure 1, at two fixed coupHngs A = 3.0 and A = 5.0. In the "strong- 
couphng" case, A = 3.0, it can be seen that the results converge exponentially 
fast in L, whereas in the "weak-coupling" regime, A — 5.0, the convergence 
is more like at these lattice sizes. This behaviour merits some further 

explanation. 

A similar phenomenon occurs in the case of the U(l) theory in (2-|-1)DE1IeI 
In the strong-coupling regime, where the mass gap is large, the usual exponen- 
tial convergence occurs. In the weak-coupling regime, however, where the mass 
gap M is very small, the finite-size scaling behaviour for small lattice sizes is 
that of a massless theory, and it is only at much larger lattice sizes L sa 1/M 
that a crossover to exponential convergence occurs. An "effective Lagrangian" 
corresponding to free, massless gluons (non-interacting QCD) should describe 
the finite-size behaviour in the present case, in line with the idea of asymptotic 
freedom. By analogy with the (2-|-l)D case, we expect a dependence for 
the corrections to the ground-state energy per site. We hope to pursue this 
analysis further at a later date. 

An anomalous feature in Figure 3b) is that the L = 8 point lies well out of 
line with the others. This occurs at other couplings also. We suspect that the 
results for L = 8 are not reliable, and that the trial wave function will have to 
be further improved to give reliable results for such large lattices. 

We have made estimates of the bulk limit, extrapolating mainly from the 
smaller L values where possible. The estimates for the bulk ground-state energy 
per site are graphed as a function^f coupling in Figure 2, where they are 
compared with previous estimates □ obtained by an 'Exact Linked Clust£r 
Expansion' (ELCE) procedure, and with the asymptotic weak-coupling seriesEj 
The Monte Carlo results agree very well with the ELCE estimates, and appear 
to match nicely onto the expected weak-coupling behaviour for A > 6. 

3.2 Wilson Loops 

The forward- walking method was used to estimate values for the m x n Wilson 
loops, W{m, n). A graph of the 'mean plaquette' 1) versus the variational 
parameter c is shown in Figure 3. A problem is immediately apparent. The 
estimate for W{1, 1) is not independent of c, in fact it depends linearly on c over 
this range, and the size of the variation is such that the probable systematic 
error due to the choice of c is an order of magnitude larger than the random 
statistical error in the results. Thus it would be advantageous in future studies 
to put more effort into improving the trial wave function, rather than merely 
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improving the statistics. 

The finite-size behaviour for the Wilson loops is similar to that of the 
ground-state energy. The estimates for the mean plaquette in the bulk limit 
are graphed as a function of coupling xLixL Figure 4, and compared with series 
estimates at strong and weak couplingaO The agreement is quite good. 



3.3 String Tension 

Having obtained estimates for the Wilson loop values on the bulk lattice, one 
can extract estimates for the 'spacelike' string tension using the Creutz ratios: 

Ka^ ~ i?„ = - In 

The results are shown in Figure 5. Also shown in Figure Sj^are some 
previous estimates derived from the 'axial' string tension, obtained lI using an 
'Exact Linked Cluster Expansion' (ELCE) method. The axial string tension aT 
is calculated as an energy per link, and must be coiHijerted to a dimensionless, 
'spacelike' tension by dividing by the 'speed of light'll3 c. We have also used the 
weak-coupling relationship between the scales of Euclidean and Hamiltonian 
lattice Yang-Mills theory calculated by Hasenfratz et alll3 to plot the results 
against the Euchdean coupling (3 — Q/g"^. 

It can be seen that the present GFMC results are in rough agreement with 
the axial string tension results in the region 4 < /? < 5, which is also the region 
where the 'roughening' transition occurs in the string tensionP For /3 > 5, 
however, the Creutz ratio i?2 runs above the ELCE estimate, and shows no 
sign of the expected crossover to an exponentially decreasing scaling behaviour 
at /3 ~ 6. We presume that this is a finite-size effect, and that the Creutz ratios 
Rn for larger n will show a substantial decrease in the 'weak-coupling' regime 
(3 > 6. That is certainly the pattern seen in the Euclidean calculations, or in 
the U(l)2+i modelQ Unfortunately, however, our present results for the larger 
Wilson loops are not of sufficient accuracy to allow worthwhile estimates of i?„ 
for n>2. 



W{n,n)W{n- 1) 
W(n,n-lV 



(13) 



4 Conclusions 

Some significant problems with the GFMC method have emerged from this 
study. The 'forward- walking' technique was introduced specifically to avoid 
any variational bias from the trial wave functionEJu As it turns out, however, 
the results for the Wilson loops show a substantial dependence on the trial 
wave function parameter c. The systematic error due to this dependence is 
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an order of magnitude larger than the statistical error, so it would pay to put 
more effort in future studies into improving the trial wave function, rather 
than simply increasing the statistics. Furthermore, the effective ensemble size 
decreases during each measurement as the descendants of each 'ancestor' state 
die out, and this produces a substantial loss in statistical accuracy at weak 
coupling, as well. 

It would be preferable if one were able to do away entirely with all the 
paraphernalia of trial wave function, weights, branching algorithms, etc, and 
just rely on some sort of Metropolis-style accept/reject algorithm to produce 
a correct distribution of walkers. Within a quantum Hamiltonian framework, 
a way is iiiown to do this, namely the Path Integral Monte Carlo (PIMC) 
approachEj We conclude that the PIMC approach may be better suited than 
GFMC to the study of large and complicated lattice Hamiltonian systems. 
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Figure 1; Ground-state energy per site graphed against 1/1/ where L is the lattice size; 
a) at coupUng A = 3.0, b) at coupling X = 5.0. The lines are merely to guide the eye. 



Figure 2: The bulk ground-state energy per site graphed against coupling A. The points 
are our Monte Carlo estimates; the solid line represents earlier ELCE estimates[19]; and the 
dashed line represents the asymptotic weak-coupling behaviour. 



Figure 3: Estimated value for the mean plaquette W{1, 1) as a function of the variational 
parameter c, for L = 6, A = 5.0. 
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Figure 4; The mean plaquette W(l, 1) for the bulk system graphed against coupling A. 
The solid line represents the strong-coupling series expansion[19], and the dashed line the 
asymptotic weaJf-coupling behaviour. 



Figure 5: The string tension Ka^ graphed against coupling /3. The circles are obtained 
from ELCE estimates of the axial string tension[19]; the triangles are Monte Carlo estimates 
of ii2. 
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